function M = slmean(X, w)
%SLMEAN Compute the mean vector or weighted mean vector(s).
%
% [ Syntax ]
%   - M = slmean(X)
%   - M = slmean(X, w)
%
% [ Arguments ]
%   - X:        The sample matrix, with each column giving a sample
%   - w:        The weights of the samples
%   - M:        The mean vector(s)
%
% [ Description ]
%   - M = slmean(X) computes the mean vector of all samples in X.
%
%     Suppose there are n samples in d-dimensional space. Then X should
%     be a d x n matrix. In the output, M is a d x 1 vector giving the
%     mean.
%
%   - M = slmean(X, w) computes the weighted mean vector(s).
%
%     w can be specified in either of the following forms:
%       - w is a 1 x n row vector, with w(i) weighting the i-th sample.
%         Then M is a d x 1 vector giving the weighted mean.
%       - w is a K x n row vector, with w(:, k) giving a group of weights
%         for all samples. Then M is a d x K matrix, with M(:, k) giving
%         the weighted mean corresponding to the weights in w(k, :).
%       - w is empty. Then it simply returns the ordinary mean vector.
%
% [ Remarks ]
%   - M should be a 2D matrix (d x n), then w should be a 1 x n row vector,
%     v would be a d x 1 column vector.
%
% [ History ]
%   - Created by Dahua Lin on Apr 22nd, 2006
%   - Modified by Dahua Lin on Sep 10th, 2006
%       - replace slmul by slmulvec to increase efficiency
%   - Modified by Dahua Lin on Jul 16th, 2007
%       - base on the new MATLAB function bsxfun
%   - Modified by Dahua Lin on Dec 19, 2007
%       - add the feature of simultaneously computing mean vectors for
%         multiple groups of weights.
%

%% parse and verify input arguments

assert(isnumeric(X) && ndims(X) == 2, ...
    'X should be a numeric matrix.');

if nargin < 2
    w = [];
else
    if isempty(w)
        k = 1;
    else
        assert(isnumeric(w) && ndims(X) == 2 && size(w,2) == size(X,2), ...
            'The weights are in invalid form.');
        k = size(w, 1);
    end    
end

%% compute

if isempty(w)     % un-weighted
    M = sum(X, 2) * (1 / size(X,2));
    
else   % weighted
    if k == 1
        w = w * (1 / sum(w));
    else
        w = bsxfun(@times, w, 1 ./ sum(w, 2));
    end
    
    M = X * w';
end


